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Equations of motion with delays naturally emerge in the analysis of complex biological control 
systems which are organized around biochemically mediated feedback interactions. We study the 
properties of a Mackey-Glass-type nonlinear map with delay - the deterministic part of the stochastic 
cerebral blood flow map (CBFM) recently introduced to elucidate the scaling properties of cerebral 
hemodynamics. We point out the existence of deterministic and nondeterministic subspaces in the 
reconstructed phase space of delay-difference equations and discuss the problem of detection of 
nonlinearities in time series generated by such dynamical systems. 
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Attractor reconstruction is certainly one of the more 
fundamental methods of analysis of nonlinear dynami- 
cal systems Q. Time series prediction, nonlinear filter- 
ing, chaos control and targeting of trajectories all exploit 
the properties of phase space Q. However, the num- 
ber of variables accessible to measurement is frequently 
limited. In many cases even the dimensionality of the 
system is not known. The fact that phase space may 
be reconstructed using the time evolution of just a sin- 
gle dynamical variable in multidimensional system paves 
the way for applications of nonlinear methods in natural 
sciences. 

As early as 1985 Babloyantz et al. demonstrated that 
certain nonlinear measures, first introduced in the con- 
text of chaotic dynamical systems, changed during slow- 
wave sleep y| . The flurry of research work that followed 
this discovery focused on the application of nonlinear dy- 
namics in quantifying brain electrical activity during dif- 
ferent mental states, sleep stages, and under the influence 
of the epileptic process (for a review see for example ) . 
Despite various technical difficulties, the number of appli- 
cations of nonlinear time series analysis has been growing 
steadily and now includes prediction of epileptic seizures 
Q , the characterization of encephalopaties Q , monitor- 
ing of anesthesia depth Q , characteristics of seizure ac- 
tivity Q , fetal ECG extraction , and control of human 
atrial fibrillation [Tcj . 

We have recently studied middle cerebral artery blood 
flow velocity in humans u sing transcranial Doppler ul- 
trasonography (TCD) [ll], We found that scaling 
properties of time series of the axial flow velocity aver- 
aged over a cardiac beat interval may be characterized by 
two exponents. The short-time scaling exponent deter- 
mines the statistical properties of fluctuations of blood 
flow velocity in short-time intervals while the Hurst ex- 
ponent describes the long-term fractal properties. To 
elucidate the nature of two scaling regions characteristic 



of healthy individuals we introduced a stochastic cere- 
bral blood flow map (CBFM). The deterministic part of 
the CBFM is reminiscent of the Mackey-Glass differential 
equation originally introduced to describe the production 
of white blood cells [T3L fl4| so that we refer to it as the 
Mackey-Glass map (MGM): 
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Herein we investigate the detailed properties of the above 
difference equation from the viewpoint of detection of 
nonlinearities and determinism. One invariably faces this 
detection problem when dealing with short, noisy physi- 
ological time series. 

It is well established that the presence of delays even in 
such dynamically simple systems as first-order differen- 
tial equations may lead to intricate and multidimensional 
evolution. Therefore, we begin our analysis with the at- 
tractor reconstruction [Til Ha, Il7| . 

Let {xi}f =l be the time series generated by Q. In the 
methods of delays the reconstructed trajectory is made 
up of the following points: 
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where J is the delay or lag and m is the embedding di- 
mension. Thus, the reconstruction depends upon two 
parameters which are not known a priori. Considerable 
effort has been devoted to finding the best approach to 
selecting judicious values of the lag J see, for example, 
(lif and references therein. This is because for finite data 
sets, especially those contaminated by noise, the choice of 
delay determines the outcome of reconstruction. In Fig. 
^we compare the autocorrelation function and mutual in- 
formation of the Mackey-Glass map time series (a = 1.3, 
b = 0.7, rj = 6). These two methods are routinely used 
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to estimate the proper reconstruction lag. Typically J is 
chosen as the delay for which the autocorrelation func- 
tion drops to a certain fraction of its initial value, e.g. 
1/e 

or the first local minimum of the mutual infor- 
mation j^jj. Thus, we obtain the lag equal to 1 and 
3, respectively. It is apparent that the mutual informa- 
tion has distinct peaks at integer multiples of the delay 
r\ and therefore is particularly well suited to detection 
of dynamical delays. The exponentially decreasing am- 
plitude of the consecutive maxima is a strong indication 
of sensitive dependence on initial condition characteristic 
of chaotic dynamics so that we proceed to calculations of 
the attractor dimension and the largest Lyapunov expo- 
nent. 

Fig. [^illustrates the process of calculation of the corre- 
lation integral C m {r, J) for embedding spaces of increas- 
ing dimension m with J = 1 |2 1| . The slope of the linear 
scaling regions in a double logarithmic plot converges to 
the correlation dimension D2 — 4.85. We obtain the 
same value for J = 3. 

In Fig. |3 we graph the average divergence of initially 
nearby trajectories as a function of the number of itera- 
tions of the MGM [H. In the scmilogarithmic graph the 
slope of the dashed straight line yields the value of the 
largest Lyapunov exponent Ai = 0.069 which confirms 
that the underlying dynamics is indeed chaotic. One can 
see in this plot that the average time after which neigh- 
boring trajectories reach the boundaries of the available 
phase space (as indicated by the constant value of their 
separation) roughly coincides with the delay time interval 
where the maxima of the mutual information are clearly 
pronounced cf. Fig. ^ 

Difference equations with delays such as the Mackey- 
Glass map are seldom discussed in standard nonlinear 
dynamics textbooks. However, let us revisit one of the 
three most fundamental nonlinear models - the Henon 
map: 
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From the second equation in @ we have y n -\ — Xn-2 
and using this expression we arrive at: 
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Hence, the two-dimensional Henon map is equivalent to 
the single difference equation Q with delay equal to 2. 
In the same vein, we can rewrite as a set of r\ ordinary 
difference equations: 
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which may be interpreted as a perturbed periodic shift 
map: 
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This linear map yields periodic orbits of length 77 regard- 
less of the initial conditions. For brevity's sake, the evo- 
lution of a difference equation with delay may be written 
with the help of the shift operator Df. 

x n+1 = D /X „ = [/(4"),4 1 )),4 1 ),...,4"- 1 )] T . (7) 

Now we turn the problem around and ask about the 
prospects of detecting determinism in short, possibly cor- 
rupted with environmental noise, time scries generated 
by a difference equation with delay, such as the MGM QJ. 
Distinguishing chaos from environmental noise has been a 
longstanding challenge to the life sciences where observed 
time series usually have relatively few data points. Sugi- 
hara and May |22j pointed out that for chaotic dynamical 
systems the accuracy of the nonlinear forecasting expo- 
nentially falls off with increasing prediction-time interval 
at a rate related to the value of the Lyapunov exponent. 
On the other hand, for uncorrelated noise the forecasting 
accuracy is roughly independent of prediction interval. 
Their original idea of using short-terrnprediction to de- 
tect determinism was later extended 153, 0, Hf| . In 
Fig. 01 we display the normalized prediction error as a 
function of the prediction time. The simplest zeroth or- 
der phase-space algorithm was employed to forecast 
the values of the MGM time series. The forecasting accu- 
racy not only initially does improve with the increasing 
prediction-time interval but also exhibits strong oscilla- 
tions. For some relatively short forecasting steps, the ac- 
curacy for the MGM time series is essentially the same as 
that of the corresponding surrogate data (data which as 
closely as possible share both the spectral properties and 
distribution with the original time series, but from which 
all nonlinear correlations have been excised [26l E^ ] ~) . To 
shed some light on this effect, let us assume that we fore- 
cast one step ahead (s = 1) and embed the MGM time 
series with 77 = 6 in three dimensions (m = 3) with delay 
J = 2 (for such delay we simply choose every other ele- 
ment of the original time series to form three dimensional 
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delay vectors (0)). For clarity of presentation we disre- 
gard the influence of the linear term in the shift operator 
Df, i.e. f(xn,Xn^) w g(xn )■ Using the projection 
operator PyfcX„ = [xn\xn , Xn ] T , we may describe the 
prediction in the reconstructed phase-space as: 

-Pl35 x ri+1 = Pl35(D g X n ) 

= P 135 [ 9 (4 6) ),4VL 2 \...,ai 5) ] T (8) 
= [g(x^),x^,x^] T = ^(P 6 24x„). 

The propagator T\ depends on variables which have been 
discarded in the embedding process so that it should not 
come as a surprise that accurate forcasting in this case is 
precluded since we selected a nondeterministic subspace 
with respect to the prediction-time interval. However, if 
we choose to forecast two steps ahead (s — 2) then 

-Pl35 x n+2 = -Pl35-Dg x n+1 

= PiM\g{x^),g{x^),x£\...,x^] T (9) 
= [g(x^),x^,x^f = T 2 {P lA ) 

and the chosen subspace is closed with respect to the 
propagator We refer to such subspace as determin- 
istic. In Fig. [3] we verify our analysis by plotting the 
normalized prediction error as a function of the neighbor- 
hood size used to construct a locally linear approximation 
of the MGM dynamics [25I Eij ] . It is apparent that for 
s = 1 the forecasting accuracy is hardly different from 
that of the corresponding surrogate data. On the other 
hand, for s = 2 the error curve with a clearly pronounced 
minimum at small neighborhood sizes is characteristic of 
nonlinear dynamical systems. The conclusions we have 
drawn are valid also for small to moderate linear coupling 
we so far ignored by setting f{x^\ Xn ) ~ g{xn^). How- 
ever, strong couplings may give rise to even more subtle 
effects which we shall discuss elsewhere. 

The forecasting error is one of many statistics applica- 
ble to testing for nonlinearities in time series. For exam- 
ple, higher-order autocorrelations such as 

< (X n - X n _ d f > J < (X n - X n ^d) 2 > (10) 

may be used to measure time asymmetry, a strong signa- 
ture of nonlinearity 0- Such autocorrelations are also 
helpful in pinpointing another interesting property of 
delay-difference equations. Upon inspection of |(BJ, one 
immediately realizes that the evolution of the perodic 
shift map amounts merely to "shuffling" of the compo- 
nents of the state vector. However, random shuffle of the 
data is the crudest way of generation of surrogate data. 
If the dimensionality of the system determined by the 
delay n is sufficiently high, it is not surprising that this 
type of time reversal test may not detect nonlinearities. 



We generated the MGM time series (a = 1.3, b = 0.7, 
w = 20) along with 39 surrogate data sets and calculated 
the corresponding value of (|TU|l . For 11 out of 20 delays d 
taken from the interval [1, 20], we could not reject the hy- 
pothesis that at the 95% confidence level the MGM time 
series results from a Gaussian linear stochastic process. 

Delay-difference equations naturally emerge in the 
analysis of complex biological control systems [11] . The 
unique properties of such dynamical systems may lead to 
difficulties in recognizing nonlinearities in their evolution 
especially when experimental time series are short and 
and delays are long. To draw attention to this issue we 
gave this paper a somewhat provocative title Determin- 
istic Uncertainty. 
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FIG. 1: Autocorrelation function (dashed line) and mutual 
information (solid line) for the time series generated by the 
MGM Q with the following parameters: a = 1.3, 6 = 0.7, 
rj = 6. 
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FIG. 2: Correlation integral Cm for embedding spaces of in- 
creasing dimension m calculated for the MGM time series 
used in Fig. 
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FIG. 3: Average divergence e of initially close trajectories as 
a function of number of iterations n of the MGM (a = 1.3, 
b — 0.7, r\ — 6) is plotted in semilogarithmic scale. The 
slope of the straight dashed line yields the value of the largest 
Lyapunov exponent Ai. The reconstruction was done with 
m = 7 and J = 1. 
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FIG. 4: The normalized forecasting error as a function of 
prediction time. The simplest zeroth order phase-space algo- 
rithm was employed to predict the values of the MGM time 
series and the corresponding surrogate data. The embedding 
was done in three dimensions with the lag J = 1. The fore- 
casting error was normalized by the standard deviation of the 
time series. 
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FIG. 5: The normalized prediction error as a function of the 
neighborhood size used to construct a locally linear approx- 
imation of the MGM dynamics (a — 1.3, b = 0.7, r\ = 6) 
or that of the surrogate data. The embedding was done in 
three dimensions with the lag 3 = 2. The displayed curves 
correspond to the length s of the prediction-time interval 



